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Abstract 

Electric field fluctuations play a major role in dissociation reactions in liquid water and determine 
its vibrational spectroscopic response. Here, we study the statistics of electric fields in liquid water 
using molecular dynamics computer simulations with a particular focus on the strong but rare 
fields that drive dissociation. Our simulations indicate that the important contributions to the 
electric field acting on OH bonds stem from water molecules less than 7 A away. Long-ranged 
contributions play a minor role. 
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I. INTRODUCTION 



The chemistry of water plays a crucial role in a wide range of processes in nature and 
technology ranging from enzymatic processes and the metabolism of living cells to acid-base 
reactions and proton transfer in hydrogen fuel cells. The theoretical study of these reactions 
is complicated by the fact that often water molecules act as reactants and products, as well 
as solvent. The archetypal example of such a complex reaction is autoionisation in liquid 
water. Here, a water molecule dissociates forming a pair of hydronium and hydroxide ions 
that then separate determining the pH of water. 

The autoionisation reaction in liquid water has been studied in experiments from which 
the reaction rate constants for dissociation and recombination have been determined 
While these experiments provided valuable information on the kinetics of the reaction, their 
space and time resolution was not sufficient to clarify the atomistic details of the reaction. In 
principle, this information can be extracted from molecular dynamics simulations with forces 
calculated ab intio via density functional theory. Such simulations are, however, complicated 
by the widely disparate time scales present in the problem; while the typical vibration 
frequency of OH bonds is of the order of tens of femtoseconds, it takes an intact water 
molecule about 10 hours to dissociate on average. This time scale problem has been recently 
overcome iSjl using a combination of Car-Parrinello molecular dynamics 141 and transition 

nn 

path sampling |5|, |6[. This approach, which does not require the definition of a reaction 
coordinate as in previous ab initio simulations of autoionisation 7|,l8l], permits to concentrate 
on the reactive event and generate unbiased dissociation trajectories. Analysis of such 
trajectories yielded the following reaction mechanism. The dissociation of the water molecule 
is initiated by a rare but strong electric field fiuctuation generated by the surrounding 
water molecules. The nascent hydronium and hydroxide ions then separate along a chain 
of hydrogen bonds via Grotthuss transfer events [9i]. When the hydrogen bond chain still 
connecting the two ions breaks, rapid recombination is prevented and from this state the 
ions can move apart further and separate completely. 

While these transition path sampling simulations produced a detailed and physically 
appealing picture of the dissociation process, they did not identify the origin of the electric 
field fiuctuation involved in the first stages of the dissociation process. In particular, no 
local hydrogen bonding pattern or ion coordination number could account for the presence 
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of the bond-breaking electric field fluctuation. Also, it remained uncertain whether these 
electric fields, sufficiently strong to initiate the dissociation, stem from water molecules near 
the dissociating one or if long ranged effects are important. 

In this paper we study electric field fluctuations in liquid water using the TIP4P model 
[lo| and investigate whether long-range electrostatics or specific local configurations are 
responsible for the strong but rare electric fields involved in the autoionisation process. In 
this particular model, water molecules are treated as rigid and the electrostatic interactions 
between molecules are modeled by point charges placed on the molecules. Since the model 
does not describe the cleavage and formation of covalent bonds, the dissociation event cannot 
be studied directly. Nevertheless, hydrogen bonds and long ranged electrostatic interactions 
are reproduced correctly such that information on the statistics of electric fields can be 
gleaned from our simulations. 

Recently, electric field fluctuations have also been considered in the context of the vi- 



brational spectroscopy of liquid water UJ, [l^, lisl, [l^. In such experiments, vibrational 



frequencies, sensitive to the local molecular environment, are used as probes for the dynam- 
ics of the hydrogen bond network percolating through the liquid. Using molecular dynamics 
simulations combined with a perturbative treatment of a simplified model of the OH oscil- 
lator, Geissler and collaborators [12I have shown that the frequency of the OH-stretch of a 
HOD molecule in liquid D2O almost perfectly correlates with the electric field acting on the 
hydrogen site. An analysis of the microscopic molecular configurations revealed that the 
strength of the electric field, and hence the frequency of the OH-stretch, is mainly deter- 
mined locally by the geometry of the HOD molecule with respect to the acceptor molecule; 
other molecules from the first solvation shell play only a minor role. At short times, vibra- 
tional decoherence is dominated by changes in this local geometry and only for longer times 
do collective density and polarization fluctuations, unrelated to specific molecular motions, 
determine the spectroscopic response. 

The remainder of this article is organised as follows. In Sec. [TTl we described the simu- 
lation methodology and explain how we calculate and classify electric fields acting on OH 
bonds. Results are presented and discussed in Sec. Ullland conclusions are given in Sec. [IVl 
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II. METHODS 



A. Simulations 

As we are interested in electric field contributions from remote hydration shells, we need 
to study systems with at least ^ 10^ water molecules. Currently, ab initio simulations with 
systems sizes of this magnitude are unfeasible, such that we have to resort to empirical 
water models. While such models do not describe bond breaking and forming, they capture 
effects due to long-ranged electrostatic interactions. In our simulations, we used the TIP4P 



interaction potential 



10|, which reproduces the structure and dynamics of liquid water 



reasonably well, even compared to more recent potentials [l5|. The very good agreement 
between experiment and simulation obtained in Ref. [isl for the spectroscopic response of 
HOD in liquid D2O, which is closely related to the statistics of electric fields, confirms the 
suitability of empirical potentials to model electric field fluctuations in liquid water. 

The TIP4P model consists of four intermolecular interaction sites placed at fixed rela- 
tive positions. Two positive point charges are placed at the hydrogen atom centres and a 
compensating negative charge is placed on an extra point on the bisector of the HOH angle. 
In addition, a Lennard- Jones (LJ) interaction site is placed at the oxygen site. The total 
intermolecular interaction potential for two TIP4P water molecules i and j with respective 
charges and and whose oxygen atoms are separated by the distance Rij, can be written 
as 

v^^ = Ty^^+4^--^^ (1) 

a p t:'Ic.O I iaj/3 J^ij J^ij 

where Viajp is the distance between the charge a on molecule i and the charge {3 on molecule 
j. The parameters of the TIP4P model are given in Tab. [H 

All data presented in this paper were obtained from a 5 ns molecular dynamics simulations 
of 1000 water molecules performed at temperature 300 K in a cubic box with side length 
L ^ 31.04 A corresponding to a density p = 1.0 g/cm^. The equations of motion were 



integrated using the velocity Verlet 16| scheme with a time step of 0.5 fs and the RATTLE 
algorithm [l^ was used to constrain bond lengths and keep the molecules rigid. The long 



ranged electrostatic interactions were calculated with Ewald summation 
space cut-off of half the box size. 



18l | with a real 
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Symbol Parameter 


Value 


Q, 


HOH angle 


104.52 


roH 


oxygen-hydrogen distance 


0.9572 A 


roM 


oxygen-extra point distance 


0.15 A 


QO 


extra point charge 


-1.04 e 


m 


hydrogen charge 


+0.52 e 


A 


repulsive LJ parameter 


0.6 A12 kcal/mol 


C 


attractive LJ parameter 


610.0 A^ kcal/mol 



TABLE I: Parameters of the TIP4P potential. 



B. Electric field calculation 

As we are interested in the electric fields involved in the dissociation of a water molecule, 
we study the electric fields at the centre of the molecule's OH bonds. Other researchers 



have studied e^ 
hydrogen site 



ectric field fluctuations at the center of charge of the molecule Uj and at the 



isl . The latter choice is particularly relevant for the dynamics of OH-stretch, 
because the force acting on the OH-bond is a mass weighted sum of the forces on the oxygen 
atom and the hydrogen atom. Due to the smaller mass of the hydrogen atom, this sum 
is dominated by the hydrogen force, motivating the choice to study electric fields at the 
hydrogen site. We have, however, verified that the electric fields at the midpoint of the OH- 
bond and at the hydrogen position, while different in magnitude, are strongly correlated. 
In particular, our simulations indicate that unusually large electric fields occurring in these 
positions originate in the same molecular configurations. 

We place Cartesian coordinates in such a way, that the x-axes lie in direction of the 
OH bonds (see Fig. [1]). A large positive value of the electric field in x-direction, E^, will 
thus correspond to a force that pulls the hydrogen atom away from the oxygen atom. For a 
particular configuration, the electric field on the OH bond is calculated from the electrostatic 
forces acting on test charges located at the centre of the OH bonds, obtained from Ewald 
summation. The distributions of the electric field components are presented in Fig. [2l Only 
Ez, the component perpendicular to the molecular plane, has a vanishing average, while 
both Ex and Ey have a positive average value. The field distributions are Gaussian to 
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FIG. 1: Cartesian coordinates located on the OH bonds. The origins he on the centres of the OH 
bonds (roH/2 ~ 0.4786 A). The basis vectors he in direction of the OH bond, pointing towards 
the respective hydrogen atom. The basis vectors Uy he also in the molecular plane, whereas the 
basis vectors are perpendicular to the plane. (Visualisations of water molecules and simulation 
snapshots presented in this article were made with the VMD software package {3].) 



a high degree, indicating that the fields result from a sum of many contributions. The 
averages (E-i) and the standard deviations (j{Ei) for the distribution functions of the three 
field components E^, Ey, and Ez are presented in Tab. [Til It is interesting to note that 
the strong electric fields observed in the simulation generate forces of similar magnitude as 
the average forces acting during the dissociation event as computed by Trout and Parrinello 
using ab initio simulations [3, S]- A direct comparison, however, is not possible since the 
rigid empirical potential used here does not permit the calculation of forces during the 
dissociation event, which is accompanied by a reorganisation of the charge distribution. 




-0.04 



-0.02 



0.02 
E. [a.u. ] 



0.04 



0.06 



FIG. 2: Distributions p{Ei) of the electric field components E^, Ey, and Ez at the centre of the 
OH bonds (see Fig. I]), 1 a.u. k. 51.422 V/A. 



The temporal persistence of electric field fiuctuations can be inferred from the normalised 
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time autocorrelation functions 

_ (E.(O)E.(t)) - {E,y 
^'^'^ ~ {ED - {E.y 

depicted in Fig. [31 Here, the angular brackets denote equilibrium averages. The time 

autocorrelation functions of the three electric field components show distinct behaviour. In z- 

direction, perpendicular to the molecular plane, the autocorrelation function has pronounced 

oscillations with a period of about 40 fs, probably due to librational motions 20|. The 



x-component shows a slightly longer decorrelation time than the other two components, 
but after ^ 250 fs the autocorrelation functions of all three electric field components decay 
exponentially with a time constant r ^ 500 fs. Our results are in good agreement with those 
of Mukamel and collaborators, who calculated electric field distributions and autocorrelation 
functions at the centre of charge of the molecule for different interaction potentials including 
TIP4P 



X 


y 


z 


(Ei) 0.0237 


0.0166 


0.0000 


a{Ei) 6.77x10-3 


6.18x10-3 


7.76x10-3 



TABLE II: Gaussian fit parameters in atomic units for the distribution functions of the electric 
field components E^, Ey, and Ez on the OH bonds (see Fig. [2]). 



C. Shell Contributions 

In this section, we study how the observed electric field distributions are generated by the 
charges of the surrounding water molecules. To this end we introduce the shell contribution 
concept. 

The total electric field Ej^ on test site a of molecule i receives contributions from all 
charges qjjs from molecules j 7^ i at distances Viajf^, 

= 7T— H H 4^riaj/3. (3) 

The component of Ejq, along the OH bond is 

E = E, = ^,^- n^^, (4) 

7 
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FIG. 3: Normalised autocorrelation functions of the electric field components at the centre of 
the OH bonds (see Fig. [1]). In the inset, the same autocorrelation functions are plotted with a 
logarithmic scale on the y-axis. 



where nf^ denotes the unit vector in the direction of the OH bond at test site a on molecule 
i. In the following we consider only the field component in x-direction and we omit the 
subscript x for simplicity. 

We now consider concentric spherical shells of width Ar surrounding the test site on 
the OH bond (see Fig. H]). Every charge qi that is closer to the test site than the cut-off 
distance, will be located in one of these shells. For a given system configuration, some 
shells will contain one or more charges, whereas others will be empty. The electric field 
contribution Es{j) of a certain shell j to the x-component of the field at the test site is 
calculated as the sum over all electric field contributions Ei from charges qi located in this 
shell, 

i?3(j) = Ei?.(r.)A,(r,). (5) 

1=1 

Here, is the position of charge qi and A'^c the total number of charges within the cut-off 
distance. The characteristic function Aj(r) is used to determine if a charge at position r is 
located within shell j, 

1 1 if r is in shell j, 
^^■(^) = io else. 

In this way, the total electric field E can be written as a sum over all shell contributions. 



E = Y.Esij), 



(7) 
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where A^s the total number of shells. 

We note that while the time evolution of the system is computed with electrostatic 
forces calculated via Ewald summation, direct summation of the electric fields is used for 
determining the shell contributions. For the system size considered here, the electrostatic 
forces obtained via Ewald and real space summation differ only slightly (see Tab. IIIip . 




FIG. 4: Shell contributions: in this two-dimensional illustration, the electric field on the centre 
of an OH bond receives contributions from positive (+) and negative (~) charges of the solvating 
water molecules. Es{j) is the electric field originating from charges located in a spherical shell j. 
The shell j is the volume (shaded region) between two concentric spheres with radii rj and rj + Ar. 

To study the origin of the large electric field fluctuations, we calculate the distributions 
Pj{Es\E) of the electric field contribution Eg of shell j, given that the observed total electric 
field E has a certain value E, 

P (E \E) = ms-E,U,r^)]6[E-E{r^)]) 

^ {S[E-E{r^)]) ■ ^ ^ 

Here, denotes a configuration of the system including the positions of all atoms. From 

this distribution, we then evaluate the average contribution Es{j,E) of shell j, for a given 

observed electric field E, 

{E^{j)S[E-Eir^)]) 



Esim 



{6[E-E{r^)]) 

dE,P,{E,\E)E,. (9) 



In the computer simulation, rather than imposing a sharp value E for the total electric field 
E, we determine the distributions for values of E that are in a narrow interval {E, E + AE). 

parameter real space Ewald sum relative difference 

(E) 0.0239 0.0237 1% 

a{E) 7.48x10-3 6.77x10-^ 10% 

TABLE 111: Comparison of the Gaussian fit parameters (in atomic units) for the distribution 
functions p{E), obtained from shell contributions in real space and from the Ewald sum for a 
system size of 1000 molecules. 

III. RESULTS AND DISCUSSION 

The average shell contributions for different values of the total field are depicted in Fig. [51 
In this calculation we considered A^^s = 200 concentric shells, with width Ar ^ 0.075 A. The 
largest shell radius corresponds to half the box size, L/2 ^ 15.52 A. Rather than using the 
shell index j = 1, . . . ,Ns, we denote the shells by their radial distance r from the test site. 
The peaks in the curves of Fig. [5] become more pronounced and are shifted to smaller radii 
for larger total fields. For all considered values of the total electric field E, however, the 
non-zero contributions to the field come only from shells less than 7 A away from the test site 
on the OH bond centre. These correspond to the contributions of charges from molecules 
belonging to the first and second hydration shell of the water molecule hosting the test site. 
Hydration shells further away have no effective net contribution to the total field on the OH 
bond. This observation remains true even for large positive values of E, suggesting that rare 
electric fields involved in autoionisation in liquid water are a rather local effect. 

The distributions of shell contributions, Pj{Es\E), for various values of the total electric 
field are shown as colour coded maps in Fig. [6] for weak (a), average (b), and strong (c) 
electric fields. For very small observed values of the total field E [Fig. El^a)], the average 
shell contributions are close to zero, even for small shell distances. Typical values of the 
total field [Fig. [6](b)] are the result of positive contributions from shells at 1.7 to 2.7 A, 
and negative contributions from shells at 2.4 to 3.3 A. Shells further away yield only a 
small positive net contribution. Larger values of E [Fig. [6](c)] are generated in a similar 
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£5" 

-0.01 
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FIG. 5: Average contributions Es{r\E) from charges in shells at distance r to the total electric 
field E at the OH bond centre, plotted as functions of r for different ranges of the total electric 
field E (values in parenthesis next to the line colour key). 

way, but with larger and more strongly peaked positive and negative shell contributions. 
The contribution peaks are shifted towards smaller shell distances, yielding higher absolute 
values of E. For all values of E, shells at distances larger than 7 A show no significant 
contributions. 

A. Configuration analysis 

Sharp peaks in the distributions Pr{Es\E), best visible in the case of large positive values 
of the total electric field E [Fig. [6](c)], indicate that typical arrangements exist where one or 
more charges at a certain distance produce a specific shell contribution Eg. For small dis- 
tances, the attribution of large electric fields to particular configurations of water molecules 
is geometrically possible. As the averaged contributions Es{r\E) rapidly decrease to values 
close to zero for shell distances r > 4 A (Fig. [5]), in our analysis we will only consider 
molecules with at least one charge closer than 4 A to the test site of the central molecule 
(molecule A in Fig. [7j). 

In most configurations that yield strong total electric fields in the range 0.050 a.u < E < 
0.055 a.u. the hydrogen atom next to the test site is involved in a hydrogen bond to 
another water molecule (B in Fig. [7]). The sharp peaks for shell radii between 2 and 3 A 
are caused by charges on this water molecule accepting the hydrogen bond. The peak at 
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distance 2.1 - 2.4 A and fields of 0.05 - 0.07 a.u. [region I in Fig. EJ^c)] originates from the 
negative charge go- The other peak at distance 2.3 - 3.0 A and fields of approximately 
—0.02 a.u. [region II in Fig. EJ^c)] is due to the charge qn of the hydrogen atom on the same 
molecule. Note that both of the curved ridged shapes created by the peaks of adjacent 
shells follow the 1/r^ dependence of the electric field. The predominance of molecule B in 
determining electric field fluctuations acting on the OH-bond was already pointed out in Ref. 
[13I . As a consequence, vibrational spectroscopy of the OH-stretch is relatively insensitive 
to the motion of water molecules other than B and thus provides only a rather local probe 
of hydrogen bond dynamics. 

Many configurations have one or two molecules (Di and D2) donating hydrogen bonds to 
the oxygen atom of the molecule (A) hosting the test site. The diffuse peak at distance 1.8 
- 2.2 A and fields of 0.01 - 0.03 a.u. [region III in Fig. [6t^c)] can be attributed to the charges 
gn on the hydrogen atoms of these molecules. 

The charges on a water molecule C accepting a hydrogen bond from the hydrogen atom 
opposite the test site show no significant contribution to the electric field E. 

The peaks in the distributions Pr{Es\E) for strong fields, caused by the charges of the 
hydrogen bonded water molecule B adjacent to the test point, are essentially responsible for 
the positive and negative peaks at shell distances 2.2 - 2.5 A and 2.5 - 3.0 A, respectively, 
in the averaged contributions Es{r\E) depicted in Fig. [51 The shoulder in the positive peak 
stems from the hydrogen charges of water molecules (Di and D2) donating a hydrogen bond 
to A. 

B. Local hydrogen bond patterns 

Since the shell contributions can be traced back to the different types B, C, and D 
of molecules hydrogen bonded to the molecule A carrying the test site, we calculated the 
probabilities of the different local hydrogen bonding patterns for configurations yielding high 
values of the electric field and compared them to those configurations with average and low 
values of the electric field (see Tab. IIV|) . We denote the local hydrogen bond pattern of a 
molecule A by the number of molecules B, C and D bonded to it. For instance, '110' denotes 
a configuration with one hydrogen bond from A to B, one hydrogen bond from A to C, and 
no bonds from molecules Di and D2 to A. Other hydrogen bond patterns are designated 
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analogously. The geometric criterion chosen to detect hydrogen bonds in this analysis is 
that the oxygen-oxygen distance roo is less than 3.5 A and the angle between ron of 
the donor and roo is less than 30°. The specific figures in Tab. IIVI depend on this choice 



quantitatively, but not qualitatively. Using these criteria, the average number of 



lydrogen 



bonds per water molecule is 3.52 as reported in earlier studies of TIP4P water 2l[ |. 



Strong fields involve a higher number of hydrogen bonds. In this case, the molecule on 
average engages in 4.27 hydrogen bonds, with 1.79 bonds donated and 2.49 accepted. This 
refiects in patterns 112, 113, and 102 occurring most frequently. Configurations with less 
than three hydrogen bonds are not detected here. Average fields involve 3.54 hydrogen 
bonds, with 1.77 equally accepted and donated. Here, the most frequent configurations are 
112 and 111. Except 112, patterns involving four or more hydrogen bonds are less frequent 
than in case of strong fields. In contrast, configurations with three hydrogen bonds or less 
become more frequent. Weak electric fields involve on average 2.38 hydrogen bonds, with 
1.32 accepted and 1.06 donated. In this case, the configurations 111 and Oil are encountered 
most often. Patterns involving four or more hydrogen bonds occur only rarely, whereas 
configurations with two or less hydrogen bonds become frequent. An analysis carried out 
for equilibrium configurations of liquid water [l4|, albeit with a different water model and a 
different hydrogen bond criterion, yielded similar frequencies of particular hydrogen bonding 
patterns as our analysis in the case of average electric field strengths. 

It is interesting to note that for strong electric fields, in a large fraction of all configurations 
(more than 40%) the water molecule with the test site accepts 3 hydrogen bonds. To a 
smaller degree, this over-coordination occurs also for typical fields, where about 5% of all 
water molecules accept more than two hydrogen bonds. While the relative weight of such 
configurations depends on the particular hydrogen bond criterion used in the analysis, our 
results clearly indicate that anomalous hydrogen bonding patterns play an important role 
in the generation of large electric fields. Whether this behavior is a refiection of the real 
situation or is due to a deficiency of the empirical force field used in our simulation is 
an interesting question that can be answered computationally only using more accurate 
potential energy surfaces such as those obtained from electronic structure calculations. 

Going from strong to weak fields, the average number of hydrogen bonds per molecule 
drops significantly from 4.27 to 2.38. In addition, the ratio of donated to accepted hydrogen 
bonds shifts: for strong fields, a molecule on average donates less hydrogen bonds than it 
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accepts, whereas for weak fields the opposite can be observed. 

In summary, large electric fields on the test site strongly correlate with a higher than 
average number of hydrogen bonds and more bonds being accepted than donated by the 
water molecule carrying the test site. 



configuration 




probability for 




BCD 


weak field 


average field strong field 





0.013 


n non 


n non 


1 


0.045 


0.005 


0.000 


2 


0.012 


0.010 


0.001 


1 


0.074 


0.000 


0.000 


1 1 


0.293 


0.013 


0.000 


1 2 


0.110 


0.056 


0.000 


1 


0.016 


0.003 


0.000 


1 1 


0.033 


0.051 


0.002 


1 2 


0.010 


0.075 


0.130 


1 3 


0.000 


0.006 


0.076 


1 1 


0.073 


0.005 


0.000 


1 1 1 


0.222 


0.210 


0.003 


1 1 2 


0.081 


0.499 


0.394 


1 1 3 


0.004 


0.051 


0.367 


1 1 4 


0.000 


0.001 


0.019 



TABLE IV: Probabilities of local hydrogen bond patterns, for configurations yielding weak (0.000 
- 0.005 a.u.), average (0.022 - 0.026 a.u.), and strong (0.050 - 0.055 a.u.) electric fields along the 
OH bond. The local configurations are defined by the number of different types of water molecules 
B, C and D hydrogen bonded to the molecule A on which the electric field is measured (see Fig. [7]). 
Patterns with relative frequencies of less than 1% are not shown. 
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IV. CONCLUSION 



We have studied the distribution functions of the electric field, generated by solvent water 
molecules, on the centres of the OH bonds of a water molecule. In order to understand the 
origin of electric field fluctuations that are sufficiently strong to initiate the dissociation of 
a water molecule, we expressed the value of the electric field component along the bond 
as the sum of electric field contributions from charges in spherical shells around the centre 
of the OH bonds. By calculating the average shell contributions for different values of the 
observed field component, we were able to conclude that the relevant contributions come 
from shells at distances smaller than 7 A, corresponding to charges on water molecules of 
the first and second hydration shell, even for the largest values of the electric field observed 
in the simulation. This suggests that the mechanism leading to the autoionisation of water 
does not involve long-ranged rearrangements in the liquid. Indeed, particularly strong fields 
correlate with an increase in the number of hydrogen bonds to the molecule's neighbours 
and a higher number of hydrogen bonds being accepted than donated by the molecule on 
which the OH bond is considered. 

Since the electric field acting on the water molecule has been shown to be almost uniquely 
related to the frequency of the OH-stretch pj^] , exceptionally strong fields may be accessible 
to experimental probes. Analysis of the respective molecular arrangements may then yield 
information on the likelihood of molecular patterns that drive chemical reactions through 
strong electric field fluctuations. 

All simulations discussed in this paper have been performed for rigid water molecules in- 
teracting via an empirical potential. The fields acting on a water molecule that is permitted 
to respond to the perturbation and actually undergo a dissociation are likely to be different 
from those reported here. This may be the reason why local hydrogen bond patterns are 
found to correlate with strong fields in our simulations, but not in Ref. [3|. Further simu- 
lations with dissociable water molecules will be necessary to clarify this issue. Such studies 
require an ab initio calculation of forces and are becoming feasible on current computer 
systems for moderate system sizes. 
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FIG. 6: Shell contribution distributions Pr{Es\E) as a function of shell distance r and field contri- 
bution Eg for (a) weak total electric fields (0.000 - 0.005 a.u.), (b) typical total electric fields (0.020 
- 0.025 a.u.) and (c) strong total electric fields (0.050 - 0.055 a.u.). The peaks of ^^.(^s = 0\E) 
for small shell distances r, caused by frequent absence of charges in shells with small volume, 
are omitted. The peaks in regions I - III in graph (c) correspond to contributions from specific 
arrangements of water molecules as explained in detail in the main text. 
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FIG. 7: A typical configuration of water molecules generating a strong electric field, 0.050 a.u. 
< E < 0.055 a.u., at the centre of the OH bond of the central molecule (dot on molecule A). 
Hydrogen bonds between water molecules are indicated by the dotted red lines. Molecule B accepts 
a hydrogen bond from molecule A involving the hydrogen next the test site. Molecule C accepts 
a hydrogen bond via the other hydrogen of molecule A. Molecules Di and D2 donate hydrogen 
bonds to molecule A. 
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